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Abstract 

Ab initio molecular dynamics is used to study deformations of sodium clus- 
ters at temperatures 500- • • 1100 K. Open-shell Nai4 cluster has two shape 
isomers, prolate and oblate, in the liquid state. The deformation is stabilized 
by opening a gap at the Fermi level. The closed-shell Nas remains magic also 
at the liquid state. 
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Electronic shell structure of alkali metal clusters has been established in numerous ex- 
periments (for a review see Ref. This shell structure corresponds to single-electron 
levels in a spherical potential and can be understood in terms of the simple jellium model 
0. It is also well-known that clusters with an open electronic shell are deformed from the 
spherical shape. In theory, the deformation has been extensively studied with jellium-type 
models f||-|5| , and it is also the result of ab initio calculations and simple tight-binding 
models |]9|,[T0|] . These calculations seem to indicate that the geometries of smallest clusters 
(N < 20) are determined solely by the symmetries of the single electron wavefunctions and 
are insensitive to the model. Moreover, even the small nuclei with equal number of nucleons 
seem to have the same geometries as alkali metal clusters fllT,5] . 



The experimental information of the deformation is obtained from the intensity varia- 
tion of the mass spectrum, or from the photoabsorption experiments []l] where the observed 
energies and the intensity ratios of the plasmon peaks are interpreted in terms of defor- 
mations of the valence electron distribution [p~2] , p~3|] . In these experiments the clusters are 
known to be hot, very likely liquidlike ]H , |I5|] . While the temperature can be incorporated 



in the electronic degrees of freedom of the jellium model []nj, so far there has not been 
attempts to systematically study the interplay between the deformation of the ionic system 
and electronic structure at elevated temperatures by detailed dynamical simulations. 

In this paper we report results obtained from ab initio molecular dynamics (MD) simula- 
tions of hot (500 < T < 1100 K ) "magic" Na 8 and open-shell Na M clusters. Our calculations 
confirm the magic nature of Nag by showing a large HOMO-LUMO gap in the proposed T=0 
ground state geometry (dodecahedron) JFj. This gap remains clear also at T rj 550 K, where 
the cluster is already liquidlike. In the open-shell Nan the deformation opens a gap, which 
also remains in the liquid state. By studying the cartesian moments of inertia, calculated 
from the coordinates of ions, we show that Na§ remains the most spherical of our clusters 
at high temperatures. In contrast, liquidlike Nai4 cluster favour axially deformed (prolate 
and oblate) shapes, and the deformation of the ionic system is driven by the opening of 
the HOMO-LUMO gap, which is directly demonstrated by starting a simulation from a 



practically spherical Nai4 isomer. 



Since the ab initio MD method we use is fully documented elsewhere [F7|, we only sum- 
marize here the key ingredients of our calculations. We perform the electronic-structure 
calculations in the framework of the density-functional theory, with the local density ap- 
proximation for the exchange-correlation energy. The valence electron - ion interaction is 
described by norm-conserving pseudopotentials |19] with a plane-wave basis. For a given 



ionic configuration the Kohn-Sham one-electron equations are solved using an iterative ma- 
trix diagonalization scheme. From the converged solution, the Hellman-Feynman forces on 
the ions are calculated, and these forces together with the classical Coulomb forces be- 
tween the positive ion cores determine the ion trajectories. We remark that due to the 
self-consistent solution of the electronic structure for each ionic configuration, the electrons 
are kept strictly on the Born-Oppenheimer potential energy surface, consequently, the length 
of the time step in the molecular dynamics simulation is limited only by the vibrational time 
scale of ions and the performance of the integration algorithm, just as in ordinary classical 
MD. In this study we use a fifth-order predictor-corrector algorithm with a time step of 5.2 
fs, giving a good conservation of the (generalized) constant of motion in the MD simulation 
[TEfl. We also remark that the method does not employ periodic boundary conditions to the 



ionic system, i.e., no supercells are involved in the calculation. For that reason the method is 



particularly suitable for studies of finite charged or multipolar systems fL? |. The calculations 
are performed typically on a cubic cell with a dimension of 45 a.u. and using a plane-wave 
kinetic energy cutoff E c = 4.4 Ry. Test calculations using larger cutoffs up to 20.1 Ry show 
that by using Ec = 4.4 Ry the ionization potential of a sodium atom deviates 1 % from the 
converged value of 5.21 eV (experimental value 5.139 eV [plj]), and the dissociation energy 
and bond length of sodium dimer are within 2 % and 5 %, respectively, from the converged 
values 0.82 eV and 2.98 A(experimental values 0.80 eV and 3.08 A, respectively f2"Tf ). 

The initial configurations for the MD simulation are prepared by setting the ions to the 
desired geometry with a reasonable nearest-neighbour distance (~ 3.5 A) and optimizing the 
structure to a local minimum on the potential energy surface. We wish to emphasize here 



that for our purposes it is not necessary to find an absolute global ground state, rather we 
wish to start the simulations from cluster isomers with representative electronic properties 
and (in the case of Nai4) desired deformations in the shape of the ionic configuration. For 
Nas, we consider only the dodecahedron. It has been found to be the ground state in previous 
pseudopotential-LDA calculations, where the ionic and electronic structure are optimized 
according to the Car-Parrinello scheme f7j. Our optimized dodecahedron has electronic 
properties very similar to those found in Ref. @, showing a large gap of 1.1 eV between 
the highest occupied and the lowest unoccupied one-electron level. The angular momentum 
character of the four occupied levels is easily identifiable showing the lowest level to be s- 
type and the next three levels p-type, of which two higher ones are degenerate. The lower 
p-level is split from them by 0.4 eV and the gap between the s- level and the lowest p- level 
is 1.2 eV. Although a different structure (stellated tetrahedron) was found in all-electron 
configuration-interaction calculations || , it has been shown |7j that all the low-lying isomers 
of Na§ have similar electronic properties, and the cluster geometry we have selected serves 
well as an example of a magic cluster having a closed electronic shell. By looking at the 
principal moments of inertia, we find that two of them are degenerate and have a larger 
value than the third one, thus the cluster can be characterized as prolate, the ratio of the 
two values I z /I r (I r degenerate, I z non-degenerate moments) being 0.76. 

We consider three isomers for Nai4, the shape of which can be characterized as prolate, 
oblate and "spherical". The prolate isomer has the lowest energy. Its geometry was chosen 
according to a recent study using Hiickel method ||10|| . The cluster consists of three stacked 
squares, the middle one rotated 45° with respect to the others. The squares are capped by 
an atom at both ends of the long axis. The ratio I z /I r is 0.39. Of the seven occupied one- 
electron states, two lowest (identified as s and p) and two highest (d-type, degenerate) states 
show a single angular momentum component, the three middle ones having some degree of 
p-d mixing. The HOMO-LUMO gap is clear, 0.4 eV. The gap is greatly reduced for the two 
other isomers which results in a fractional occupancy around the Fermi level (see Ref. [|T8|]). 
The geometric structure of the oblate isomer is obtained by relaxing two stacked hexagonal 

4 



layers of seven atoms (one atom in the middle), rotated 30° with respect to each other. It is 
interesting to note that the energy of this rather artificial structure is only 0.3 eV (or 0.02 
eV/atom) higher than that of the prolate isomer. It is probable that we could find oblate 
isomers which are energetically even closer to the prolate isomer, paralleling the results of 
Lauritsch et al. Q, who find oblate and prolate shapes (of the electron density distribution) 
to be separated only by 0.1 eV in their jellium calculations for Na^. Finally, the "spherical" 
isomer, having degenerate moments, is obtained by relaxing a structure where an atom is 
placed on one (100) face of the 13-atom cuboctahedron. The energy of the relaxed cluster 
is 0.9 eV higher than that of the prolate isomer. Looking at the electronic structure we find 
that the Fermi level lies in the middle of split ld-levels, mixed with the 2s (according to the 
nomenclature in the spherical potential). 

Next we follow the evolution of the geometric and electronic structure of each cluster 
in molecular dynamics runs, where the cluster is given enough kinetic energy to raise the 
(vibrational) temperature |22| to 500-600 K, after which the system is allowed to evolve 
at constant energy. The rapid heating results in a complete disorder (melting) during the 
first 1-2 ps, which is easily verified by monitoring the mean-squared displacements (msd) of 
atoms as a function of time. After the melting, a typical linear rise in msd versus time is 
observed, and the diffusion constant estimated from the slope of the msd curve falls in the 
range (2 — 9) x 10 -5 cm 2 /s, not far from the values (9 — 13) x 10 -5 cm 2 /s expected for bulk 
liquid sodium at these temperatures |23[ . 



The fluctuations of the geometric shape of the clusters are shown in Fig. 1, which shows 
the evolution of the three cartesian moments of inertia as a function of time, and in a 
two-dimensional "shape space", where a particular time step tk is represented by a point 
{h{tk)/h(tk), h(tk)/ h(tk)), where I\ < I 2 < I-s- Trajectories in this space are confined 
within an upper-left triangle of a square, with corners of (0,0), (0,1), and (1,1). All the 
ideal prolate shapes (ii < I 2 = h) fall on the line (0,1) — ^(1,1) while the ideal oblate shapes 
(ii = I 2 < I3) form the diagonal (0,0)— >(1,1) of the square. The corner (1,1) obviously has 
all the shapes with degenerate moments, i.e. shapes which have at least the cubic symmetry. 



As can be seen from Fig. la, the trajectory for Nas probes relatively uniformly the region 
around a triaxial average shape, characterized by coordinate (0.72,0.89). The evolution of 
the electronic structure is shown in Fig. 2a, where we plot the Kohn-Sham eigenenergies 
as a function of time. We see only minor changes compared to zero-temperature electronic 
structure: the HOMO-LUMO gap averages to 0.9 eV, and the splitting of the p-levels 
increases to 0.5 eV. It is interesting to note the survival of the zero-temperature gap also 
between the (unoccupied) ld-2s manifold and the If levels (only two of them are plotted). 
This indicates that the average triaxial shape of the ionic configuration does not destroy too 
much the (approximately) spherical background potential felt by the electrons. The overall 
behaviour of the electronic shell structure confirms the magic nature of Nag even at T=550 
K, in accordance with previous Car-Parrinello calculations [[/]]. 

We turn now to the more interesting case of the open-shell Nai4 cluster (see figs, lb-d 
and 2b-d). The results shown for the principal moments of inertia (Fig. lb,c) for the runs 
starting from the prolate and oblate shapes show that both of these shapes are indeed stable 
(meaning that there is a marked energy barrier between them). Furthermore, transitions 
between isomers appear to be possible in the time scale of our simulations, t « 10 ps. (Fig. 
lb). The gap visible at T=0 for the prolate isomer persists also in the finite-temperature 
simulation (Fig. 2b), and there is a gap opening for the oblate isomer as soon as the dynamics 
is turned on (Fig. 2c). It is instructive to note that even the Kohn-Sham eigenvalues reflect 
the shape of the cluster: for prolate shapes, the lowest p-type orbital (corresponding to 
the state along the longer axis) is separated from the higher p-orbitals whence for oblate 
shapes there is a gap between the two lowest (states along the two long axes) and the higher 
p-orbital. In both cases we observe a considerable mixing between the higher p-orbital(s) 
and the d-orbitals. Finally, Fig. 2d shows a "dynamical Jahn- Teller effect" taking place in 
the "spherical" isomer of Nai4. Within 0.5 ps after the dynamics is turned on, we see the 
splitting of the ld-2s manifold and the gap opening at ep. The shape of the cluster deforms 
to the prolate side (Fig. Id), hence also the lowest p-state splits off from the two higher 
p-states. 
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By passing we want to mention that when the Na i4 cluster is heated to 1100 K the 
time-averaged level density does not any more show a pronounced gap at the Fermi level. 
At this high temperature cluster shape changes continuously between all shapes and the 
trajectory in the "shape space" covers uniformly a large area. It should be noted that a free 
cluster, initially at 1100 K, would very fast cool by evaporation, though evaporation is not 
observed in our simulations due to the short time scale. Thus it is likely that at experimental 
temperatures the Nai4 cluster is strongly deformed also in the liquid state. 

In conclusion, we have compared the behaviour of hot " magic" Na 8 and open-shell Nai 4 
clusters in ab initio molecular dynamics simulations at elevated vibrational temperatures 
from 500 to 1100 K, where the ionic structure of all the clusters is completely disordered 
and exhibits liquidlike diffusion. Na 8 remains magic around 500 K showing a clear electronic 
shell structure. We find two stable shape isomers (prolate and oblate) for Na i4 and observe 
inter-isomeric transitions in our 10 ps time scale for clusters around 600 K. We demonstrate 
a "dynamical Jahn- Teller effect" in a spherical isomer of Nai 4 where the Fermi level initially 
lies in the middle of the ld-2s manifold, but once the dynamics is turned on, the cluster 
opens a gap at ep, which stabilizes the deformation to the prolate shape. 
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FIGURES 

FIG. 1. Time evolution of the cartesian moments of inertia of liquid sodium clusters, (a) Nas 
at 550 K, (b) prolate Nai4 at 650 K, (c) oblate Nai4 at 610 K, (d) Nai4 at 480 K started from a 
spherical shape. For each cluster we show the trajectory in the phase space and the three moments 
of inertia (in units of 10 3 (amu)ag) as a function of the time in ps (inset). Note the transition from 
the prolate shape to oblate shape (and back) in (b), and the rapid strong deformation in (d). 

FIG. 2. Evolution of single-electron eigenvalues (in eV) as a function of time in ps (left) and 
the time-averaged density of states (in arbitrary units, right) in (a) Nag at 550 K, (b) prolate Nai4 
at 650 K, (c) oblate Nai4 at 610 K, (d) Nai4 at 480 K started from a spherical shape. Note the 
gap opening at ep in (d). 



11 



